Lista 2

Questão 1

O seu professor chega na sala de aula e mostra uma moeda. Você suspeita que a moeda possa ser falsa e ter duas caras. Considere a priori probabilidades iguais para os eventos da moeda ser falsa ou ser honesta (i.e. uma moeda bem equilibrada).

Calcule a sua probabilidade de obter cara num lançamento dessa moeda.

Para calcular a probabilidade de obter “cara” em um lançamento dessa moeda, considerando que a moeda pode ser honesta ou falsa. Suponha que:

$H$ seja o evento “a moeda é honesta”.
$F$ seja o evento “a moeda é falsa”.
$C_{i}$ seja o evento “obter cara no i-ésimo lançamento”.

Assumindo pelo enunciado:

A priori, as probabilidades de a moeda ser honesta ou falsa são iguais: $P(H) = P(F) = 0.5 $.
Se a moeda é honesta, a probabilidade de obter cara é $P(C_{1}|H) = 0.5$.
Se a moeda é falsa, a probabilidade de obter cara é $P(C_{1}|F) = 1$ (pois tem duas caras).

A probabilidade de obter cara será dada por:

$P(C_{1}) = P(C_{1}|H)P(H) + P(C_{1}|F)P(F)$

Substituindo os valores fornecidos:

$P(C_{1}) = (0.5 \times 0.5) + (1 \times 0.5)$

$P(C_{1}) = 0.25 + 0.5$

$P(C_{1}) = 0.75$

Portanto, a probabilidade de obter cara em um lançamento dessa moeda é $0.75$ ou 75%.

Se o professor lançar a moeda e o resultado for cara, qual é agora a probabilidade dela ser falsa?

Queremos encontrar $P(F|C_{1})$, a probabilidade de a moeda ser falsa dado que o resultado foi cara. Segue via Teorema de Bayes:

$P(F|C_{1}) = \frac{P(C_{1}|F)P(F)}{P(C_{1})}$

Logo:

$P(F|C_{1}) = \frac{P(C_{1}|F)P(F)}{P(C_{1})} = \frac{1 \times 0.5}{0.75} = \frac{0.5}{0.75} = \frac{2}{3}$

Portanto, a probabilidade de a moeda ser falsa, dado que o resultado do lançamento foi cara, é $\frac{2}{3}$ ou aproximadamente 66.67%.

Se o professor lançar a moeda n vezes e obter n caras, qual é a probabilidade dela ser falsa? Estude o comportamento desta probabilidade para n grande.

Assumimos:

$C^n$ o evento “obter cara em todos os $n$ lançamentos”.

Queremos calcular $P(F|C^n)$, a probabilidade de a moeda ser falsa dado que todos os $n$ lançamentos resultaram em cara.

Segue via Teorema de Bayes:

$P(F|C^n) = \frac{P(C^n|F)P(F)}{P(C^n)}$

Onde:

$P(F) = 0.5$ (a priori, a moeda é falsa).
$P(H) = 0.5$ (a priori, a moeda é honesta).
$P(C^n|F) = 1$ (se a moeda é falsa, a probabilidade de obter cara em todos os lançamentos é 1).
$P(C^n|H) = 0.5^n$ (se a moeda é honesta, a probabilidade de obter cara em todos os $n$ lançamentos é $0.5^n$).

Para encontrar $P(C^n)$, recorremos a probabilidade total:

$P(C^n) = P(C^n|H)P(H) + P(C^n|F)P(F)$

Substituindo os valores fornecidos:

$P(C^n) = (0.5^n \times 0.5) + (1 \times 0.5)$

$P(C^n) = 0.5^{n+1} + 0.5$

Agora, substituímos na fórmula de Bayes:

$P(F|C^n) = \frac{1 \times 0.5}{0.5^{n+1} + 0.5}$

$P(F|C^n) = # Inferência Bayesiana {.unnumbered}

Observa-se que valores grandes de $n$, $0.5^n$ se aproxima de 0 muito rapidamente, porque $0.5^n$ decresce exponencialmente.

Portanto, para $n$ grande, $P(F|C^n)$ se aproxima de:

$P(F|C^n) \approx \frac{1}{0 + 1} = 1$

Isso significa que, conforme $n$ aumenta, a probabilidade de a moeda ser falsa se aproxima de 1 (ou 100%).

Portanto, se o professor lançar a moeda $n$ vezes e obter $n$ caras, a probabilidade de a moeda ser falsa tende a 1 conforme $n$ se torna grande.

Se o professor lançar a moeda uma vez e obter cara, qual é a probabilidade do pr´oximo lançamento ser cara?

Após um lançamento resultando em “cara”, a probabilidade de a moeda ser falsa, $P(F|C)$, foi calculada como $\frac{2}{3}$. Consequentemente, a probabilidade de a moeda ser honesta, $P(H|C)$, será:

$P(H|C) = 1 - P(F|C) = 1 - \frac{2}{3} = \frac{1}{3}$

Para a probabilidade de obter “cara” no próximo lançamento, dado que o primeiro lançamento foi “cara” basta obter a soma das probabilidades condicionais de obter “cara” no próximo lançamento, dadas as duas possibilidades sobre a moeda (falsa ou honesta):

$P(C_2|C_1) = P(C_2|F)P(F|C_1) + P(C_2|H)P(H|C_1)$

Onde:

$P(C_2|C_1)$ é a probabilidade de obter “cara” no próximo lançamento, dado que o primeiro lançamento foi “cara”.
$P(C_2|F) = 1$ (se a moeda é falsa, a probabilidade de obter “cara” no próximo lançamento é 1).
$P(C_2|H) = 0.5$ (se a moeda é honesta, a probabilidade de obter “cara” no próximo lançamento é 0.5).
$P(F|C_1) = \frac{2}{3}$ (a probabilidade de a moeda ser falsa após o primeiro “cara”).
$P(H|C_1) = \frac{1}{3}$ (a probabilidade de a moeda ser honesta após o primeiro “cara”).

Substituindo os valores, temos:

$P(C_2|C_1) = (1 \times \frac{2}{3}) + (0.5 \times \frac{1}{3})$

$P(C_2|C_1) = \frac{2}{3} + \frac{1}{6}$

$P(C_2|C_1) = \frac{4}{6} + \frac{1}{6}$

$P(C_2|C_1) = \frac{5}{6}$

Portanto, a probabilidade de obter “cara” no próximo lançamento, dado que o primeiro lançamento resultou em “cara”, é $\frac{5}{6}$ ou aproximadamente 83.33%.

Explique porque é falso neste contexto a afirmação ”os dois lançamentos da moeda são independentes”, e explique qual seria a afirmação correta.

No contexto descrito, a afirmação “os dois lançamentos da moeda são independentes” é falsa porque o resultado de cada lançamento influencia nossa crença sobre a natureza da moeda (se ela é honesta ou falsa). Essa influência altera a probabilidade dos resultados subsequentes.

Vale ressaltar que dois eventos $A$ e $B$ são independentes se e somente se:

$P(A \cap B) = P(A) \cdot P(B)$

Ou, de forma equivalente:

$P(B|A) = P(B)$

Se os lançamentos fossem verdadeiramente independentes, o resultado de um lançamento não alteraria a probabilidade do resultado do próximo lançamento. No entanto, neste contexto, o resultado do primeiro lançamento afeta nossa crença sobre a natureza da moeda, que por sua vez afeta a probabilidade do segundo lançamento. Ou seja , os lançamentos são dependentes porque $P(C_2 | C_1) \neq P(C_2)$. Aqui, $C_1$ e $C_2$ representam o primeiro e o segundo lançamento resultando em cara, respectivamente.

A afirmação correta seria: “Os dois lançamentos da moeda são condicionalmente dependentes”.

Em resumo, o resultado de um lançamento influencia a probabilidade do próximo lançamento devido à incerteza sobre a natureza da moeda. Isso demonstra que os lançamentos não são independentes; são, na verdade, condicionalmente dependentes.

Questão 2

Seja $y_{1} , y_{2} ,\dots, y_{n}$ uma amostra da distribuição de Bernoulli com probabilidade de sucesso θ e considere uma distribuição a priori uniforme para θ.

Ache a distribuição a posteriori de θ e a sua média e variância.

Para uma amostra de Bernoulli, a verossimilhança é dada por:

$P(y_1, y_2, \ldots, y_n | \theta) = \theta^{\sum_{i=1}^n y_i} (1 - \theta)^{n - \sum_{i=1}^n y_i}$

Seja $k = \sum_{i=1}^n y_i$, que é o número de sucessos na amostra. Assim, a função de verossimilhança pode ser escrita como:

$P(y_1, y_2, \ldots, y_n | \theta) = \theta^k (1 - \theta)^{n - k}$

A distribuição a priori uniforme para $\theta$ no intervalo $[0, 1]$ é dada por:

$P(\theta) = 1 \quad \text{para} \quad \theta \in [0, 1]$

Distribuição a posteriori:

$P(\theta | y_1, y_2, \ldots, y_n) \propto P(y_1, y_2, \ldots, y_n | \theta) \cdot P(\theta)$

Substituindo as funções:

$P(\theta | y_1, y_2, \ldots, y_n) \propto \theta^k (1 - \theta)^{n - k} \cdot 1$

$P(\theta | y_1, y_2, \ldots, y_n) \propto \theta^k (1 - \theta)^{n - k}$

Essa é a forma da distribuição Beta. Portanto, a distribuição a posteriori de $\theta$ é uma distribuição Beta com parâmetros $\alpha = k + 1$ e $\beta = n - k + 1$:

$\theta | y_1, y_2, \ldots, y_n \sim \text{Beta}(k + 1, n - k + 1)$

Para uma distribuição $\text{Beta}(\alpha, \beta)$, a média e a variância são dadas por:

$\mathbb{E}[\theta] = \frac{\alpha}{\alpha + \beta}$

$\text{Var}(\theta) = \frac{\alpha \beta}{(\alpha + \beta)^2 (\alpha + \beta + 1)}$

Substituindo os parâmetros $\alpha = k + 1$ e $\beta = n - k + 1$:

Média:

$\mathbb{E}[\theta | y_1, y_2, \ldots, y_n] = \frac{k + 1}{(k + 1) + (n - k + 1)} = \frac{k + 1}{n + 2}$

Variância:

$\text{Var}(\theta | y_1, y_2, \ldots, y_n) = \frac{(k + 1)(n - k + 1)}{((k + 1) + (n - k + 1))^2 ((k + 1) + (n - k + 1) + 1)}$

$\text{Var}(\theta | y_1, y_2, \ldots, y_n) = \frac{(k + 1)(n - k + 1)}{(n + 2)^2 (n + 3)}$

Mostre que é possı́vel expressar a esperança a posteriori de θ da forma $(1−w)\times E(θ)+w\times \hat{θ̂}$, onde E(θ) e θ̂ são respectivamente a esperança a priori e a estimativa máximo verossı́mil de θ, e interprete este resultado.

A distribuição a posteriori de $\theta$ foi encontrada como $\text{Beta}(k + 1, n - k + 1)$, tendo $\mathbb{E}[\theta | y_1, y_2, \ldots, y_n] = \frac{k + 1}{n + 2}$.

A priori, temos uma distribuição uniforme para $\theta$ no intervalo $[0, 1]$, cuja esperança é:

$\mathbb{E}[\theta] = \frac{1}{2}$

Sabendo que a estimativa de máxima verossimilhança (MLE) de $\theta$ para uma amostra Bernoulli é a proporção de sucessos na amostra:

$\hat{\theta} = \frac{k}{n}$

Basta substituir os valores encontrados:

$\mathbb{E}[\theta | y_1, y_2, \ldots, y_n] = \frac{k + 1}{n + 2}$

$\mathbb{E}[\theta] = \frac{1}{2}$

$\hat{\theta} = \frac{k}{n}$

Vamos expressar $\frac{k + 1}{n + 2}$ na forma desejada:

$\frac{k + 1}{n + 2} = \left(1 - \frac{2}{n + 2}\right) \times \frac{1}{2} + \frac{2}{n + 2} \times \frac{k}{n}$

Para verificar, expandimos a expressão:

$\left(1 - \frac{2}{n + 2}\right) \times \frac{1}{2} + \frac{2}{n + 2} \times \frac{k}{n}$

$= \frac{1}{2} - \frac{1}{n + 2} + \frac{2k}{n(n + 2)}$

$\frac{k + 1}{n + 2} = \frac{n}{2(n + 2)} + \frac{2k}{n(n + 2)}$

Efetuando mais uma manipulação algébrica:

$\frac{k + 1}{n + 2} = \frac{1}{2} \left( \frac{n}{n + 2} \right) + \left( \frac{2}{n + 2} \right) \times \frac{k}{n}$

Logo, $w = \frac{2}{n + 2}$ e $(1 - w) = \frac{n}{n + 2}$. Isso significa que a esperança a posteriori de $\theta$ é uma média ponderada entre a esperança a priori e a estimativa de máxima verossimilhança, com os pesos dependendo do tamanho da amostra $n$.

Ou seja:

Quando a amostra é pequena ($n$ pequeno), a esperança a priori tem um peso maior na esperança a posteriori.
Quando a amostra é grande ($n$ grande), a estimativa de máxima verossimilhança tem um peso maior, e a esperança a posteriori se aproxima da MLE.

Isso reflete a intuição de que com mais dados, nossa crença sobre $\theta$ é mais fortemente influenciada pela evidência observada do que pelas crenças a priori.

Se $y_{n+1}$ é uma observação futura deste processo de Bernoulli, ache a distribuição preditiva $p(y_{n+1} |y_{1} , y_{2} ,\dots, y_{n})$.

Distribuição preditiva de $y_{n+1}$:

Para obter a distribuição preditiva de $y_{n+1}$, precisamos calcular a probabilidade de $y_{n+1}$ condicionado aos dados observados $y_1, y_2, \dots, y_n$. A distribuição de $y_{n+1}$ dado $\theta$ é Bernoulli com parâmetro $\theta$:

$p(y_{n+1} | \theta) = \theta^{y_{n+1}} (1 - \theta)^{1 - y_{n+1}}$

A distribuição preditiva é obtida integrando $\theta$ fora desta expressão, ponderando pela distribuição a posteriori de $\theta$, encontrada previamente:

$p(y_{n+1} | y_1, y_2, \dots, y_n) = \int_0^1 p(y_{n+1} | \theta) p(\theta | y_1, y_2, \dots, y_n) \, d\theta$

Substituindo $p(y_{n+1} | \theta)$ e $p(\theta | y_1, y_2, \dots, y_n)$:

$p(y_{n+1} = 1 | y_1, y_2, \dots, y_n) = \int_0^1 \theta \cdot \frac{\theta^k (1 - \theta)^{n - k}}{B(k + 1, n - k + 1)} \, d\theta$

$p(y_{n+1} = 0 | y_1, y_2, \dots, y_n) = \int_0^1 (1 - \theta) \cdot \frac{\theta^k (1 - \theta)^{n - k}}{B(k + 1, n - k + 1)} \, d\theta$

Para $y_{n+1} = 1$:

$p(y_{n+1} = 1 | y_1, y_2, \dots, y_n) = \int_0^1 \theta \cdot \frac{\theta^k (1 - \theta)^{n - k}}{B(k + 1, n - k + 1)} \, d\theta$

$= \frac{1}{B(k + 1, n - k + 1)} \int_0^1 \theta^{k + 1 - 1} (1 - \theta)^{n - k} \, d\theta$

$= \frac{B(k + 2, n - k + 1)}{B(k + 1, n - k + 1)}$

$= \frac{\frac{\Gamma(k + 2) \Gamma(n - k + 1)}{\Gamma(n + 3)}}{\frac{\Gamma(k + 1) \Gamma(n - k + 1)}{\Gamma(n + 2)}}$

$= \frac{(k + 1) \Gamma(k + 1) \Gamma(n - k + 1) / \Gamma(n + 3)}{\Gamma(k + 1) \Gamma(n - k + 1) / \Gamma(n + 2)}$

$= \frac{(k + 1)}{n + 2}$

Para $y_{n+1} = 0$:

$p(y_{n+1} = 0 | y_1, y_2, \dots, y_n) = \int_0^1 (1 - \theta) \cdot \frac{\theta^k (1 - \theta)^{n - k}}{B(k + 1, n - k + 1)} \, d\theta$

$= \frac{1}{B(k + 1, n - k + 1)} \int_0^1 \theta^k (1 - \theta)^{n - k + 1 - 1} \, d\theta$

$= \frac{B(k + 1, n - k + 2)}{B(k + 1, n - k + 1)}$

$= \frac{\frac{\Gamma(k + 1) \Gamma(n - k + 2)}{\Gamma(n + 3)}}{\frac{\Gamma(k + 1) \Gamma(n - k + 1)}{\Gamma(n + 2)}}$

$= \frac{\Gamma(n - k + 2) / \Gamma(n + 3)}{\Gamma(n - k + 1) / \Gamma(n + 2)}$

$= \frac{(n - k + 1)}{n + 2}$

Assim, a distribuição preditiva $p(y_{n+1} | y_1, y_2, \dots, y_n)$ é:

$p(y_{n+1} = 1 | y_1, y_2, \dots, y_n) = \frac{k + 1}{n + 2}$

$p(y_{n+1} = 0 | y_1, y_2, \dots, y_n) = \frac{n - k + 1}{n + 2}$

A distribuição preditiva de uma observação futura $y_{n+1}$ é a média ponderada das proporções de sucessos e fracassos observadas, ajustadas pelo tamanho da amostra, refletindo o princípio de suavização de Laplace, onde mesmo com uma grande quantidade de dados, há um pequeno ajuste baseado na crença a priori.

Questão 4

No exercício 2, calcule:

a estimativa bayesiana para Perda Quadrática

Para calcular a estimativa Bayesiana sob a perda quadrática, precisamos encontrar o estimador que minimiza o valor esperado da perda quadrática, geralmente expressa como $L(\theta, d) = (\theta - d)^2$, onde $\theta$ é o parâmetro verdadeiro e $d$ é a decisão ou estimativa tomada. Logo:

$\mathbb{E}[(\theta - d)^2 | y_1, y_2, \dots, y_n]$

onde $\theta$ segue a distribuição a posteriori dado os dados observados $y_1, y_2, \dots, y_n$.

Dado que $\theta$ tem distribuição a posteriori $\text{Beta}(k + 1, n - k + 1)$, a estimativa Bayesiana que minimiza a perda quadrática é a média da distribuição a posteriori, afinal, a média de uma distribuição é o estimador de mínimo erro quadrático médio (MSE).

A média da distribuição $\text{Beta}(\alpha, \beta)$ é $\frac{\alpha}{\alpha + \beta}$. Para a nossa distribuição a posteriori $\text{Beta}(k + 1, n - k + 1)$:

$\mathbb{E}[\theta | y_1, y_2, \dots, y_n] = \frac{k + 1}{(k + 1) + (n - k + 1)} = \frac{k + 1}{n + 2}$

Portanto, a estimativa Bayesiana de $\theta$ para a perda quadrática, dada a amostra de dados $[y_1, y_2, \dots, y_n]$, é:

$d^* = \frac{k + 1}{n + 2}$

A estimativa Bayesiana $d^*$ corresponde à média a posteriori do parâmetro $\theta$, levando em conta tanto a informação a priori (que $\theta$ é uniformemente distribuído entre 0 e 1) quanto a informação proveniente da amostra ($k$ sucessos em $n$ tentativas). O ajuste de $+1$ no numerador e $+2$ no denominador é uma forma de regularização que evita extremos baseados em amostras pequenas, fornecendo uma suavização estatística que protege contra conclusões extremas com dados limitados, refletindo um compromisso entre o conhecimento a priori sobre a distribuição de $\theta$ e a evidência obtida da amostra, alinhando-se com os princípios fundamentais da inferência Bayesiana.

o limite da estimativa bayesiana para Perda Zero-Um quando $\epsilon→ 0$. No caso especial que $n = 12$, $s = \sum ^{12}_{i=1} y_{i} = 9$

A função de perda zero-um é dada por $L(\theta, d) = 1$ se $\theta \neq d$ e $L(\theta, d) = 0$ se $\theta = d$. Essencialmente, essa perda penaliza decisões incorretas sem considerar quão “erradas” elas estão, ao contrário da perda quadrática que penaliza de acordo com o quadrado da diferença entre a verdadeira $\theta$ e a estimativa $d$.

Cabe ressaltar que na estimativa para perda zero-um, procuramos minimizar o erro de classificação. Em geral, o estimador que minimiza a esperança da perda zero-um sob uma distribuição de probabilidades é o valor de $\theta$ que maximiza a distribuição a posteriori, ou seja, a moda da distribuição a posteriori.

No nosso caso especial, a distribuição a posteriori de $\theta$ após observar os dados é a $\text{Beta}(s + 1, n - s + 1)$, onde $s = \sum_{i=1}^n y_i$. Para $n = 12$ e $s = 9$, temos:

$\theta | y_1, y_2, \dots, y_{12} \sim \text{Beta}(9 + 1, 12 - 9 + 1) = \text{Beta}(10, 4)$

A moda de uma distribuição $\text{Beta}(\alpha, \beta)$, para $\alpha > 1$ e $\beta > 1$, é dada por:

$\text{Mode}[\theta] = \frac{\alpha - 1}{\alpha + \beta - 2}$

Para a nossa distribuição $\text{Beta}(10, 4)$:

$\text{Mode}[\theta] = \frac{10 - 1}{10 + 4 - 2} = \frac{9}{12} = 0.75$

O limite de $\epsilon \rightarrow 0$ na pergunta sugere considerar um caso extremo da função de perda ou dos parâmetros envolvidos. Neste contexto, normalmente analisamos o comportamento da estimativa à medida que o tamanho da amostra cresce ou como a estimativa se comporta sob condições extremas da amostra ou dos parâmetros, todavia, como já calculamos a moda, o que seria o resultado para a perda zero-um, não há dependência direta de $\epsilon$ no cálculo, a menos que $\epsilon$ esteja relacionado a algum tipo de regularização ou parâmetro de suavização não especificado explicitamente na questão.

a estimativa bayesiana sob Perda Absoluta

Para encontrar a estimativa Bayesiana sob a perda absoluta, precisamos minimizar o valor esperado da perda absoluta. A função de perda absoluta é dada por $L(\theta, d) = |\theta - d|$, onde $\theta$ é o parâmetro verdadeiro e $d$ é a decisão ou estimativa tomada.

O estimador que minimiza a perda absoluta é a mediana da distribuição a posteriori, afinal, a mediana minimiza a soma das distâncias absolutas em relação a ela mesma.

A distribuição a posteriori de $\theta$ após observar os dados $ y_1, y_2, , y_n $ é $\text{Beta}(k + 1, n - k + 1)$, onde $k = \sum_{i=1}^n y_i$. No caso especial onde $ n = 12 $ e $ k = s = 9 $, temos:

$\theta | y_1, y_2, \ldots, y_{12} \sim \text{Beta}(10, 4)$

Para uma distribuição $\text{Beta}(\alpha, \beta)$, a mediana não tem uma fórmula fechada simples, mas pode ser aproximada ou calculada numericamente.

Uma aproximação comum para a mediana de uma distribuição $\text{Beta}(\alpha, \beta)$ é:

$\text{Mediana} \approx \frac{\alpha - \frac{1}{3}}{\alpha + \beta - \frac{2}{3}}$

Substituindo $\alpha = 10$ e $\beta = 4$:

$\text{Mediana} \approx \frac{10 - \frac{1}{3}}{10 + 4 - \frac{2}{3}} = \frac{9.6667}{13.3333} \approx 0.725$

Portanto, a estimativa Bayesiana de $\theta$ sob a perda absoluta, dada a amostra específica, é aproximadamente 0.725.

Ao optar por uma aproximação via python, obtemos:

Show/hide code

# Parâmetros da distribuição Beta
a, b = 10, 4

median = round(beta.median(a, b), 3)



print(f"A mediana da distribuição Beta({a}, {b}): {median}")

A mediana da distribuição Beta(10, 4): 0.725

De modo a encontrar o mesmo resultado para a estimativa Bayesiana de $\theta$ sob a perda absoluta.

um intervalo HPD com nível 99%.

O intervalo de alta densidade posterior (HPD - Highest Posterior Density) é um intervalo que contém uma determinada percentagem da densidade de probabilidade e onde todas as probabilidades dentro do intervalo são maiores ou iguais às probabilidades fora do intervalo. Para calcular um intervalo HPD de 99% para a distribuição a posteriori de $\theta$, onde $\theta$ segue uma distribuição $\text{Beta}(10, 4)$, utilizaremos o software python, pois a distribuição Beta não permite uma solução analítica simples para o intervalo HPD.

Para $n = 12$ e $s = 9$, a distribuição a posteriori de $\theta$ é $\text{Beta}(10, 4)$. A fórmula da função de densidade para a distribuição Beta é:

$f(\theta; \alpha, \beta) = \frac{\theta^{\alpha - 1} (1 - \theta)^{\beta - 1}}{B(\alpha, \beta)}$

onde $B(\alpha, \beta)$ é a função Beta, e para nosso caso específico, $\alpha = 10$ e $\beta = 4$.

Logo:

Show/hide code

# Parâmetros da Beta
a, b = 10, 4


def hpd_interval(a, b, level=0.99):
    
    lower_bound = (1 - level) / 2
    upper_bound = 1 - lower_bound
    
    lower_quantile = beta.ppf(lower_bound, a, b)
    upper_quantile = beta.ppf(upper_bound, a, b)
    
    return lower_quantile, upper_quantile


hpd_99 = hpd_interval(a, b, 0.99)

hpd_99_formatted = (round(hpd_99[0], 3), round(hpd_99[1], 3))

print(f"Intervalo HPD de 99% para Beta({a}, {b}): {hpd_99_formatted}")

Intervalo HPD de 99% para Beta(10, 4): (0.379, 0.943)

O intervalo resultante explicita a região onde a densidade da distribuição Beta é mais alta e contém 99% da probabilidade total.

Questão 8

Suponha que $(x_{1} , x_{2} , x_{3})$ dado $p_{1} , p_{2} , p_{3}$ segue uma distribuição Multinomial com parâmetros $n$ e $(p_{1} , p_{2} , p_{3})$ , onde $p_{i} ≥ 0 \space \text{e} \space p_{1} + p_{2} + p_{3} = 1$, e que, a priori, $(p_{1} , p_{2} , p_{3})$ segue uma distribuição de Dirichlet com parâmetros $(α_{1} , α_{2} , α_{3} )$:

Ache a distribuição a posteriori de $p_{1} , p_{2} , p_{3}$ e as distribuições a posteriori marginais de $p_{i}$ (i = 1, 2, 3).

A distribuição a posteriori de $(p_1, p_2, p_3)$ é uma distribuição de Dirichlet com parâmetros atualizados.

A distribuição Dirichlet a priori é:

$(p_1, p_2, p_3) \sim \text{Dirichlet}(\alpha_1, \alpha_2, \alpha_3)$

A função de verossimilhança para a distribuição multinomial é:

$P(x_1, x_2, x_3 | p_1, p_2, p_3) = \frac{n!}{x_1! x_2! x_3!} p_1^{x_1} p_2^{x_2} p_3^{x_3}$

Afinal, a probabilidade de obter exatamente $x_1$ ocorrências do resultado 1, $x_2$ ocorrências do resultado 2, e $x_3$ ocorrências do resultado 3 é dada pelo produto das probabilidades individuais elevadas ao número de ocorrências, multiplicado pelo número de permutações possíveis.

A distribuição a posteriori será:

$P(p_1, p_2, p_3 | x_1, x_2, x_3) \propto p_1^{\alpha_1 - 1} p_2^{\alpha_2 - 1} p_3^{\alpha_3 - 1} \cdot p_1^{x_1} p_2^{x_2} p_3^{x_3}$

$\propto p_1^{\alpha_1 + x_1 - 1} p_2^{\alpha_2 + x_2 - 1} p_3^{\alpha_3 + x_3 - 1}$

Portanto, a distribuição a posteriori é uma distribuição de Dirichlet com parâmetros $(\alpha_1 + x_1, \alpha_2 + x_2, \alpha_3 + x_3)$:

$(p_1, p_2, p_3) | (x_1, x_2, x_3) \sim \text{Dirichlet}(\alpha_1 + x_1, \alpha_2 + x_2, \alpha_3 + x_3)$

As distribuições marginais de $p_i$ para uma distribuição de Dirichlet são Beta. Para $p_1$:

$p_1 | (x_1, x_2, x_3) \sim \text{Beta}(\alpha_1 + x_1, \alpha_2 + x_2 + \alpha_3 + x_3)$

Similarmente, para $p_2$:

$p_2 | (x_1, x_2, x_3) \sim \text{Beta}(\alpha_2 + x_2, \alpha_1 + x_1 + \alpha_3 + x_3)$

E para $p_3$:

$p_3 | (x_1, x_2, x_3) \sim \text{Beta}(\alpha_3 + x_3, \alpha_1 + x_1 + \alpha_2 + x_2)$

Calcule as estimativas bayesianas de $p_{i}$ e de $p_{j}$ − $p_{i}$ sob Perda Quadrática (i, j = 1, 2, 3, i 6= j).

Estimativas Bayesianas de $p_i$:

Para a distribuição Dirichlet $\text{Dirichlet}(\alpha_1, \alpha_2, \alpha_3)$, a esperança de $p_i$ é:

$\mathbb{E}[p_i | x_1, x_2, x_3] = \frac{\alpha_i + x_i}{\sum_{j=1}^3 (\alpha_j + x_j)}$

Portanto, as estimativas Bayesianas de $p_1$, $p_2$, e $p_3$ são:

$\hat{p_1} = \frac{\alpha_1 + x_1}{\alpha_1 + x_1 + \alpha_2 + x_2 + \alpha_3 + x_3}$ $\hat{p_2} = \frac{\alpha_2 + x_2}{\alpha_1 + x_1 + \alpha_2 + x_2 + \alpha_3 + x_3}$ $\hat{p_3} = \frac{\alpha_3 + x_3}{\alpha_1 + x_1 + \alpha_2 + x_2 + \alpha_3 + x_3}$

Estimativa de $p_j - p_i$:

Para calcular a estimativa Bayesiana de $p_j - p_i$ sob a perda quadrática, usaremos a linearidade da esperança. A estimativa é simplesmente a diferença das esperanças de $p_j$ e $p_i$:

$\mathbb{E}[p_j - p_i | x_1, x_2, x_3] = \mathbb{E}[p_j | x_1, x_2, x_3] - \mathbb{E}[p_i | x_1, x_2, x_3]$

Portanto:

$\mathbb{E}[p_j - p_i | x_1, x_2, x_3] = \frac{\alpha_j + x_j}{\alpha_1 + x_1 + \alpha_2 + x_2 + \alpha_3 + x_3} - \frac{\alpha_i + x_i}{\alpha_1 + x_1 + \alpha_2 + x_2 + \alpha_3 + x_3}$

Podemos assim, concatenar os resultados obtidos :

A distribuição a posteriori de $(p_1, p_2, p_3)$ é $\text{Dirichlet}(\alpha_1 + x_1, \alpha_2 + x_2, \alpha_3 + x_3)$.
As distribuições marginais de $p_i$ são $\text{Beta}(\alpha_i + x_i, \sum_{j \neq i} (\alpha_j + x_j))$.
As estimativas Bayesianas de $p_i$ são $\hat{p_i} = \frac{\alpha_i + x_i}{\sum_{j=1}^3 (\alpha_j + x_j)}$.
As estimativas Bayesianas de $p_j - p_i$ são $\hat{p_j} - \hat{p_i} = \frac{\alpha_j + x_j}{\sum_{k=1}^3 (\alpha_k + x_k)} - \frac{\alpha_i + x_i}{\sum_{k=1}^3 (\alpha_k + x_k)}$.

Questão 11

É conhecido que 25% dos pacientes de um certo grupo que sofrem de enxaqueca melhoram após duas horas de serem tratados com um placebo. Para verificar se uma droga nova é melhor que o placebo, $n = 20$ pacientes foram tratados com o placebo e verificou-se que após duas horas $s = 8$ deles relataram ter melhorado. Seja $\theta$ a probabilidade de um paciente tratado com a droga nova melhorar após duas horas.

Especifique a hipótese nula $H_0$ e a alternativa $H_1$;

Hipótese nula $H_0$: $\theta = 0.25$ (a droga nova é tão eficaz quanto o placebo)
Hipótese alternativa $H_1$: $\theta > 0.25$ (a droga nova é mais eficaz que o placebo)

Usando a distribuição a priori “não informativa” $\theta \sim \text{Uniforme}(0, 1)$, calcule as chances relativas a priori e a posteriori de $H_1$ e o correspondente Fator de Bayes.

$P(\theta > 0.25)$
$P(\theta \leq 0.25)$

A distribuição a priori uniforme $\theta \sim \text{Uniforme}(0, 1)$:

$P(\theta \leq 0.25) = 0.25$

$P(\theta > 0.25) = 1 - 0.25 = 0.75$

As chances relativas a priori:

$O(H_1) = \frac{P(\theta > 0.25)}{P(\theta \leq 0.25)} = \frac{0.75}{0.25} = 3$

Sabemos que a distribuição a posteriori é $\theta | s \sim \text{Beta}(9, 13)$.

Probabilidade a Posteriori de $ > 0.25 $:

$P(\theta > 0.25 | s) = 1 - F_{\text{Beta}}(0.25 | 9, 13)$

Consequentemente:

$P(\theta \leq 0.25 | s) = F_{\text{Beta}}(0.25 | 9, 13)$

Via python:

Show/hide code

# Parâmetros da distribuição Beta
a, b = 9, 13

posterior_prob_H1 = 1 - beta.cdf(0.25, a, b)
posterior_prob_H0 = beta.cdf(0.25, a, b)

print(f"P(θ > 0.25 | s) = {posterior_prob_H1:.3f}")
print(f"P(θ ≤ 0.25 | s) = {posterior_prob_H0:.3f}")

P(θ > 0.25 | s) = 0.944
P(θ ≤ 0.25 | s) = 0.056

Chances Relativas a Posteriori:

$O(H_1 | s) = \frac{P(\theta > 0.25 | s)}{P(\theta \leq 0.25 | s)} = \frac{0.944}{0.056} = 16.810$

Fator de Bayes:

$B_{10}(s) = \frac{O(H_1 | s)}{O(H_1)} = \frac{16.810}{3} = 5.603$

Seja $d = 1$ a decisão de rejeitar $H_0$ e $d = 0$ a de não rejeitar. Considere a função de perda de Neyman para a qual é 5 vezes mais custoso rejeitar $H_0$ quando ela é verdadeira do que não rejeitar quando ela é falsa. Calcule a decisão ótima a posteriori;

Definindo a função de perda de Neyman:

$L(d = 1, \theta \in H_0) = 5$
$L(d = 0, \theta \not\in H_0) = 5$
$L(d = 1, \theta \not\in H_0) = 0$
$L(d = 0, \theta \in H_0) = 0$

Considerando a função de perda de Neyman:

$a_0 = L(d = 1, \theta \leq 0.25) = 5L(d = 0, \theta \leq 0.25) = 5a_1$

A decisão ótima é rejeitar $H_0$ se:

$P(\theta > 0.25 | s) > \frac{a_0}{a_0 + a_1} = \frac{5a_1}{5a_1 + a_1} = \frac{5}{6} = 0.833$

Como $P(\theta > 0.25 | s) = 0.944$ é maior do que $0.833$, a decisão ótima é rejeitar $H_0$.

É razoável chamar essa distribuição a priori de “não informativa” nesse problema? Se a sua resposta for negativa, sugira uma outra distribuição a priori e refaça os cálculos anteriores.

Uma distribuição uniforme pode ser vista como “não informativa” em um certo sentido, mas pode não ser a mais adequada. Para hipóteses específicas como $H_0$ e $H_1$, uma distribuição Beta poderia ser mais informativa.

A priori uniforme favorece $H_1$ no sentido que $P(\theta > 0.25) = 3P(\theta \leq 0.25)$. Uma forma de corrigir isso é usar uma distribuição Beta com $\beta = 3\alpha$, o que dá lugar a uma priori imprópria alternativa:

$\theta \sim \text{Beta}(\alpha, 3\alpha)$

Para obter uma variância maximizada, precisamos definir $\alpha \to 0$. Então:

$P(\theta) \propto \frac{1}{\theta(1 - \theta)}$

Com a priori imprópria $P(\theta) \propto \frac{1}{\theta(1 - \theta)}$, a distribuição a posteriori é proporcional a:

$P(\theta | s) \propto P(s | \theta) P(\theta) \propto \theta^s (1 - \theta)^{n - s} \cdot \frac{1}{\theta(1 - \theta)} = \theta^{s-1} (1 - \theta)^{n-s-1}$

Isto é uma distribuição $\text{Beta}(s, n-s)$.

Para $n = 20$ e $s = 8$:

$\theta | s \sim \text{Beta}(8, 12)$

Show/hide code

alpha_prior, beta_prior = 1, 3
n = 20
s = 8

a_post = alpha_prior + s
b_post = beta_prior + n - s


posterior_prob_H1_beta = 1 - beta.cdf(0.25, a_post, b_post)
print(f"P(θ > 0.25 | s) = {posterior_prob_H1_beta:.3f}")

P(θ > 0.25 | s) = 0.904

A decisão de rejeitar $H_0$ é mantida tanto com a priori uniforme quanto com a priori Beta(1, 3).

Lista 3.1

Questão 1

Seja $x_1, x_2, \ldots, x_n$ uma amostra da distribuição de Poisson com média $\theta$, e considere a priori que $\theta$ tem uma distribuição Gama com parâmetros $\alpha$ e $\beta$ (isto é, com média $\alpha\beta^{-1}$ e variância $\alpha\beta^{-2}$).

Ache a distribuição a posteriori de $\theta$ e a sua média e variância.

Dada uma amostra $x_1, x_2, \ldots, x_n$ da distribuição de Poisson com média $\theta$, a verossimilhança é:

$P(x_1, x_2, \ldots, x_n | \theta) = \prod_{i=1}^n \frac{\theta^{x_i} e^{-\theta}}{x_i!}$

$P(x_1, x_2, \ldots, x_n | \theta) = \theta^{\sum_{i=1}^n x_i} e^{-n\theta} \prod_{i=1}^n \frac{1}{x_i!}$

Considerando que $\theta$ tem uma distribuição Gama a priori com parâmetros $\alpha$ e $\beta$:

$P(\theta) \propto \theta^{\alpha - 1} e^{-\beta \theta}$

A distribuição a posteriori é:

$P(\theta | x_1, x_2, \ldots, x_n) \propto P(x_1, x_2, \ldots, x_n | \theta) P(\theta)$

Substituindo as funções de verossimilhança e a priori:

$P(\theta | x_1, x_2, \ldots, x_n) \propto \theta^{\sum_{i=1}^n x_i} e^{-n\theta} \cdot \theta^{\alpha - 1} e^{-\beta \theta}$

$P(\theta | x_1, x_2, \ldots, x_n) \propto \theta^{(\alpha + \sum_{i=1}^n x_i) - 1} e^{-(\beta + n) \theta}$

Portanto, a distribuição a posteriori é uma distribuição Gama com parâmetros:

$\alpha^* = \alpha + \sum_{i=1}^n x_i$

$\beta^* = \beta + n$

A média e variância da distribuição a posteriori são:

$\mathbb{E}[\theta | x] = \frac{\alpha^*}{\beta^*} = \frac{\alpha + \sum_{i=1}^n x_i}{\beta + n}$

$\text{Var}[\theta | x] = \frac{\alpha^*}{(\beta^*)^2} = \frac{\alpha + \sum_{i=1}^n x_i}{(\beta + n)^2}$

Mostre que é possível expressar a esperança a posteriori de $\theta$ da forma $w\bar{x} + (1-w)\alpha\beta^{-1}$, e interprete este resultado.

Sabemos que:

$\mathbb{E}[\theta | x] = \frac{\alpha + \sum_{i=1}^n x_i}{\beta + n}$

Podemos escrever:

$\mathbb{E}[\theta | x] = \frac{\alpha}{\beta + n} + \frac{\sum_{i=1}^n x_i}{\beta + n}$

Definimos $\bar{x} = \frac{1}{n} \sum_{i=1}^n x_i$, assim:

$\mathbb{E}[\theta | x] = \frac{\alpha}{\beta + n} + \frac{n\bar{x}}{\beta + n}$

Se definirmos $w = \frac{n}{\beta + n}$, temos:

$\mathbb{E}[\theta | x] = (1 - w) \frac{\alpha}{\beta} + w\bar{x}$

Portanto, podemos escrever:

$\mathbb{E}[\theta | x] = w\bar{x} + (1-w)\alpha\beta^{-1}$

Isso mostra que a estimativa a posteriori é uma combinação ponderada entre a média a priori $\alpha\beta^{-1}$ e a média amostral $\bar{x}$.

O que acontece na parte (ii) quando $\beta$ é grande com $\alpha\beta^{-1}$ fixo? Interprete!

Quando $\beta$ é grande com $\alpha\beta^{-1}$ fixo:

$w = \frac{n}{\beta + n} \approx \frac{n}{\beta}$

Se $\beta \to \infty$, então $w \to 0$, o que significa que a esperança a posteriori $\mathbb{E}[\theta | x]$ será dominada pela média a priori $\alpha\beta^{-1}$, resultando em:

$\mathbb{E}[\theta | x] \approx (1 - 0) \alpha\beta^{-1} = \alpha\beta^{-1}$

Mostre que existe um número $c$ tal que a variância a posteriori é maior do que a variância a priori sempre que $\bar{x} > c$. Ache $c$ e interprete este resultado.

A variância a posteriori é:

$\text{Var}[\theta | x] = \frac{\alpha + \sum_{i=1}^n x_i}{(\beta + n)^2}$

A variância a priori é:

$\text{Var}[\theta] = \frac{\alpha}{\beta^2}$

Queremos encontrar $c$ tal que $\text{Var}[\theta | x] > \text{Var}[\theta]$ quando $\bar{x} > c$.

$\frac{\alpha + \sum_{i=1}^n x_i}{(\beta + n)^2} > \frac{\alpha}{\beta^2}$

Multiplicando ambos os lados por $(\beta + n)^2 \beta^2$:

$(\alpha + \sum_{i=1}^n x_i) \beta^2 > \alpha (\beta + n)^2$

$\alpha \beta^2 + \sum_{i=1}^n x_i \beta^2 > \alpha (\beta^2 + 2\beta n + n^2)$

$\sum_{i=1}^n x_i \beta^2 > \alpha (2\beta n + n^2)$

Dividindo ambos os lados por $n \beta^2$:

$\bar{x} > \frac{\alpha (2\beta n + n^2)}{n \beta^2}$

$\bar{x} > \frac{\alpha (2\beta + n)}{\beta^2}$

Portanto, o valor crítico $c$ é:

$c = \frac{\alpha (2\beta + n)}{\beta^2} = \frac{\alpha}{\beta} \times (2 + \frac{n}{\beta}) = (2 + \frac{n}{\beta}) \times \mathbb{E}[\theta]$

Questão 3

Seja $x_1, x_2, \ldots, x_n$ uma amostra da distribuição de Poisson com média $\theta$, e considere a priori que $\theta$ tem uma distribuição Gama com parâmetros $\alpha = 1$ e $\beta = 1$. Ache:

A estimativa bayesiana de $\theta$ no caso de perda quadrática.

A estimativa sob perda quadrática é imediata: $\mathbb{E}[\theta | D] = \alpha^* / \beta^* = (s + 1) / (n + 1)$.

O limite da estimativa bayesiana sob perda zero-um quando $\epsilon \to 0$.

O limite da estimativa sob perda zero-um é a moda da distribuição a posteriori, que nesse caso é $(\alpha^* - 1) / \beta^* = s / (n + 1)$.

Agora, quando $n = 10$ e $\bar{x} = 1.55$, temos que $\alpha^* = 16.50$ e $\beta^* = 11$, de forma que a estimativa sob PA é a mediana $\text{med}(\theta | D) \approx 1.470$.

Para o caso $n = 10$ e $\bar{x} = 1.55$, ache:

A estimativa bayesiana sob perda absoluta.

Temos que $\alpha^* = 16.50$ e $\beta^* = 11$, de forma que a estimativa sob PA é a mediana $\text{med}(\theta | D) \approx 1.470$.

O intervalo HPD para $\theta$ com nível 95%.

Vamos implementar os cálculos em Python para verificar os resultados obtidos anteriormente e, por fim, explicitar o intervalo HPD:

Show/hide code

n = 10
x_bar = 1.55
s = n * x_bar

alpha_post = 1 + s
beta_post = 1 + n

theta_quadratic = alpha_post / beta_post
print(f"Estimativa bayesiana sob perda quadrática: {theta_quadratic:.3f}")

theta_zero_one = (alpha_post - 1) / beta_post
print(f"Limite da estimativa sob perda zero-um: {theta_zero_one:.3f}")

theta_absolute = gamma.median(alpha_post, scale=1/beta_post)
print(f"Estimativa bayesiana sob perda absoluta: {theta_absolute:.3f}")

hpd_interval = gamma.interval(0.95, alpha_post, scale=1/beta_post)
print(f"Intervalo HPD de 95%: ({hpd_interval[0]:.3f}, {hpd_interval[1]:.3f})")

Estimativa bayesiana sob perda quadrática: 1.500
Limite da estimativa sob perda zero-um: 1.409
Estimativa bayesiana sob perda absoluta: 1.470
Intervalo HPD de 95%: (0.866, 2.306)

Questão 5

Seja $x_1, x_2, \ldots, x_n$ uma amostra da distribuição Normal com média $\mu$ e variância $\phi^{-1}$ conhecida, e considere a distribuição a priori $\mu \sim \mathcal{N}(\mu_0, \tau^{-1})$.

Ache a distribuição a posteriori de $\mu$.

Dada a amostra $x_1, x_2, \ldots, x_n$ da distribuição Normal com média $\mu$ e variância $\phi^{-1}$, a verossimilhança é:

$f(x_i | \mu) = \sqrt{\frac{\phi}{2\pi}} \exp\left( -\frac{\phi}{2}(x_i - \mu)^2 \right)$

A função de verossimilhança conjunta para a amostra é:

$L(\mu | x_1, x_2, \ldots, x_n) = \left(\frac{\phi}{2\pi}\right)^{n/2} \exp\left( -\frac{\phi}{2} \sum_{i=1}^n (x_i - \mu)^2 \right)$

A priori $\mu \sim \mathcal{N}(\mu_0, \tau^{-1})$ é dada por:

$f(\mu) = \sqrt{\frac{\tau}{2\pi}} \exp\left( -\frac{\tau}{2} (\mu - \mu_0)^2 \right)$

Onde teremos a distribuição a posteriori:

$f(\mu | x_1, x_2, \ldots, x_n) \propto L(\mu | x_1, x_2, \ldots, x_n) f(\mu)$

Substituindo as expressões:

$f(\mu | x_1, x_2, \ldots, x_n) \propto \exp\left( -\frac{\phi}{2} \sum_{i=1}^n (x_i - \mu)^2 \right) \exp\left( -\frac{\tau}{2} (\mu - \mu_0)^2 \right)$

Agrupando os termos quadráticos em $\mu$:

$\propto \exp\left( -\frac{1}{2} \left[ \phi \sum_{i=1}^n x_i^2 + \tau \mu_0^2 - 2\mu \left( \phi \sum_{i=1}^n x_i + \tau \mu_0 \right) + \mu^2 (\phi n + \tau) \right] \right)$

Completando o quadrado:

$f(\mu | x_1, x_2, \ldots, x_n) \propto \exp\left( -\frac{1}{2} \left[ (\phi n + \tau)\left(\mu - \frac{\phi \sum_{i=1}^n x_i + \tau \mu_0}{\phi n + \tau}\right)^2 \right] \right)$

Portanto, a distribuição a posteriori é uma Normal com os parâmetros atualizados:

$\mu | x_1, x_2, \ldots, x_n \sim \mathcal{N}\left( \frac{\phi \sum_{i=1}^n x_i + \tau \mu_0}{\phi n + \tau}, \frac{1}{\phi n + \tau} \right)$

Mostre que é possível expressar a esperança a posteriori de $\mu$ da forma $w\bar{x} + (1 - w)\mu_0$, e interprete este resultado.

A esperança a posteriori é:

$\mathbb{E}[\mu | x_1, x_2, \ldots, x_n] = \frac{\phi \sum_{i=1}^n x_i + \tau \mu_0}{\phi n + \tau}$

Podemos reescrever essa expressão como uma combinação ponderada entre a média amostral $\bar{x}$ e a média a priori $\mu_0$. Definimos $w = \frac{\phi n}{\phi n + \tau}$:

$\mathbb{E}[\mu | x_1, x_2, \ldots, x_n] = w \bar{x} + (1 - w) \mu_0$

onde $\bar{x} = \frac{1}{n} \sum_{i=1}^n x_i$.

Este resultado mostra que a estimativa a posteriori é uma média ponderada entre a média da amostra $\bar{x}$ e a média a priori $\mu_0$, com pesos $w$ e $1 - w$ respectivamente. O peso $w$ depende do tamanho da amostra e dos parâmetros $\phi$ e $\tau$.

Se $\bar{x}_m$ é a média de $m$ observações futuras $x_{n+1}, \ldots, x_{n+m}$, condicionalmente independentes de $x_1, \ldots, x_n$, ache a distribuição preditiva $p(\bar{x}_m | x_1, \ldots, x_n)$.

Se $\bar{x}_m$ é a média de $m$ observações futuras $x_{n+1}, \ldots, x_{n+m}$, condicionalmente independentes de $x_1, \ldots, x_n$, a distribuição preditiva é:

$p(\bar{x}_m | x_1, \ldots, x_n)$

A média de $m$ observações futuras é normalmente distribuída com a mesma média $\mu$ e variância $\frac{1}{\phi m}$. Portanto, a média preditiva é uma média ponderada da distribuição posterior de $\mu$ e da distribuição das observações futuras.

A variância total é a soma das variâncias individuais:

$\text{Var}(\bar{x}_m | x_1, \ldots, x_n) = \text{Var}(\mu | x_1, \ldots, x_n) + \text{Var}(x_{n+1}, \ldots, x_{n+m} | \mu)$

$= \frac{1}{\phi n + \tau} + \frac{1}{\phi m}$

A distribuição preditiva de $\bar{x}_m$ é:

$\bar{x}_m | x_1, \ldots, x_n \sim \mathcal{N}\left( \frac{\phi n \bar{x} + \tau \mu_0}{\phi n + \tau}, \frac{1}{\phi n + \tau} + \frac{1}{\phi m} \right)$

Discuta o que acontece com os resultados anteriores quando a distribuição a priori $p(\mu) \propto 1$ (i.e., o caso limite quando $\tau \to 0$).

Quando $\tau \to 0$, a priori se torna não informativa, ou seja, $p(\mu) \propto 1$. Nesse caso, o peso $w$ na expressão da esperança a posteriori se torna 1, pois:

$w = \frac{\phi n}{\phi n + \tau} \approx 1$

Portanto, a esperança a posteriori se torna a média amostral:

$\mathbb{E}[\mu | x_1, x_2, \ldots, x_n] \approx \bar{x}$

E a variância a posteriori se torna:

$\text{Var}(\mu | x_1, x_2, \ldots, x_n) \approx \frac{1}{\phi n}$

Isso mostra que, com uma priori não informativa, a distribuição a posteriori é dominada pelos dados da amostra, e a média a posteriori é simplesmente a média amostral.

Questão 6

(a) Ache o estimador bayesiano de $\mu$ no caso de (i) perda quadrática $(L(\hat{\mu}, \mu) = (\hat{\mu} - \mu)^2)$, (ii) perda absoluta $(L(\hat{\mu}, \mu) = |\hat{\mu} - \mu|)$ e (iii) perda zero-um $(L(\hat{\mu}, \mu) = I(|\hat{\mu} - \mu| \geq \epsilon))$.

Para uma função de perda quadrática $L(\hat{\mu}, \mu) = (\hat{\mu} - \mu)^2$, o estimador bayesiano é a média da distribuição a posteriori.

Como sabemos da questão anterior que a distribuição a posteriori é uma Normal com parâmetros atualizados:

$\mu | x_1, x_2, \ldots, x_n \sim \mathcal{N}\left( \frac{\phi \sum_{i=1}^n x_i + \tau \mu_0}{\phi n + \tau}, \frac{1}{\phi n + \tau} \right)$

A média dessa distribuição é:

$\hat{\mu}_{\text{PQ}} = \mathbb{E}[\mu | x_1, x_2, \ldots, x_n] = \frac{\phi \sum_{i=1}^n x_i + \tau \mu_0}{\phi n + \tau}$

Para a distribuição normal, a média, a mediana e a moda são todos iguais devido à sua simetria. Assim, os estimadores sob perda quadrática (i), perda absoluta (ii) e perda zero-um (iii) são todos iguais a:

$\hat{\mu} = \frac{\phi \sum_{i=1}^n x_i + \tau \mu_0}{\phi n + \tau}$

Essa equivalência é uma consequência direta da simetria da distribuição normal, onde a distribuição é perfeitamente simétrica em torno da média, que coincide com a mediana e a moda.

(b) Ache o intervalo HPD para $\mu$ com nível $100(1 - \alpha)\%$.

O intervalo HPD (Highest Posterior Density) para uma distribuição Normal $\mathcal{N}(\mu_{\text{post}}, \sigma^2_{\text{post}})$ é encontrado nos quantis da distribuição que correspondem ao nível de confiança desejado.

Os limites do intervalo HPD de $100(1 - \alpha)\%$ são:

$\left[ \mu_{\text{post}} - z_{\alpha/2} \sigma_{\text{post}}, \mu_{\text{post}} + z_{\alpha/2} \sigma_{\text{post}} \right]$

onde $\mu_{\text{post}} = \frac{\phi \sum_{i=1}^n x_i + \tau \mu_0}{\phi n + \tau}$ e $\sigma^2_{\text{post}} = \frac{1}{\phi n + \tau}$.

Vamos calcular esses valores:

Show/hide code

n = 10
x_bar = 1.55
phi = 1
tau = 1
mu_0 = 0

sum_x = n * x_bar
mu_post = (phi * sum_x + tau * mu_0) / (phi * n + tau)
sigma_post = np.sqrt(1 / (phi * n + tau))

alpha = 0.05
z = norm.ppf(1 - alpha / 2)
hpd_interval = [mu_post - z * sigma_post, mu_post + z * sigma_post]

print(f"Intervalo HPD de 95%: [{hpd_interval[0]:.3f}, {hpd_interval[1]:.3f}]")

Intervalo HPD de 95%: [0.818, 2.000]

Lista 3.2

Questão 1

Considere o modelo $X \sim \text{Uniforme} (\theta, \theta + 1)$ $(-\infty < \theta < \infty)$.

(a) Mostre que $\theta$ é um parâmetro de locação.

Um parâmetro $\theta$ é um parâmetro de locação se, ao adicionar uma constante $c$ a $\theta$, a distribuição resultante for uma translação da distribuição original.

Para o modelo $X \sim \text{Uniforme} (\theta, \theta + 1)$:

\[ f_X(x; \theta) = \begin{cases} 1 & \text{se } \theta \leq x \leq \theta + 1 \\ 0 & \text{caso contrário} \end{cases} \]

Se adicionarmos uma constante $c$ a $\theta$, a nova distribuição será $X \sim \text{Uniforme} (\theta + c, \theta + c + 1)$:

\[ f_X(x; \theta + c) = \begin{cases} 1 & \text{se } \theta + c \leq x \leq \theta + c + 1 \\ 0 & \text{caso contrário} \end{cases} \]

Essa é uma translação da distribuição original, mostrando que $\theta$ é um parâmetro de locação.

(b) Dada a priori não informativa usual para o modelo de locação, $p(\theta) \propto 1$, e uma amostra aleatória $X_1, \ldots, X_n$ do modelo acima, discuta a propriedade da distribuição a posteriori.

A priori não informativa usual para um modelo de locação é $p(\theta) \propto 1$.

Dada uma amostra aleatória $X_1, \ldots, X_n$ do modelo $\text{Uniforme} (\theta, \theta + 1)$:

A verossimilhança é:

$L(\theta | X_1, \ldots, X_n) \propto \prod_{i=1}^n f_X(x_i; \theta)$

\[ L(\theta | X_1, \ldots, X_n) = \begin{cases} 1 & \text{se } \theta \leq \min(X_i) \text{ e } \max(X_i) \leq \theta + 1 \\ 0 & \text{caso contrário} \end{cases} \]

Dado que a priori é $p(\theta) \propto 1$, a distribuição a posteriori é proporcional à verossimilhança:

$p(\theta | X_1, \ldots, X_n) \propto L(\theta | X_1, \ldots, X_n)$

\[ p(\theta | X_1, \ldots, X_n) = \begin{cases} \frac{1}{1 - (\max(X_i) - \min(X_i))} & \text{se } \max(X_i) - 1 \leq \theta \leq \min(X_i) \\ 0 & \text{caso contrário} \end{cases} \]

Portanto, a distribuição a posteriori de $\theta$ é uma distribuição uniforme no intervalo $[ \max(X_i) - 1, \min(X_i)]$.

(c) Na situação da parte (b), ache:

os estimadores bayesianos sob Perda Quadrática e sob Perda Absoluta;

Perda Quadrática: Para a função de perda quadrática $L(\hat{\theta}, \theta) = (\hat{\theta} - \theta)^2$, o estimador bayesiano é a média da distribuição a posteriori.

Dado que a distribuição a posteriori é uniforme no intervalo $[ \max(X_i) - 1, \min(X_i)]$, a média é:

$\hat{\theta}_{\text{PQ}} = \frac{(\max(X_i) - 1) + \min(X_i)}{2}$

Perda Absoluta: Para a função de perda absoluta $L(\hat{\theta}, \theta) = |\hat{\theta} - \theta|$, o estimador bayesiano é a mediana da distribuição a posteriori.

A mediana de uma distribuição uniforme no intervalo $[ \max(X_i) - 1, \min(X_i)]$ é:

$\hat{\theta}_{\text{PA}} = \frac{(\max(X_i) - 1) + \min(X_i)}{2}$

o intervalo HPD para $\theta$ com credibilidade $100(1 - \alpha)\%$;

Para a distribuição uniforme no intervalo $[ \max(X_i) - 1, \min(X_i)]$, o intervalo HPD $100(1 - \alpha)\%$ é simplesmente o intervalo central que cobre $100(1 - \alpha)\%$ da distribuição:

$\text{HPD} = \left[ X_n - 1 + \frac{\alpha}{2} (X_1 - (X_n - 1)), X_1 - \frac{\alpha}{2} (X_1 - (X_n - 1)) \right]$

Onde $\min(X_i) = X_1$ e $\max(X_i) = X_n$

a probabilidade a posteriori do evento $\theta \geq \theta_0$.

Para encontrar a probabilidade a posteriori do evento $\theta \geq \theta_0$:

$P(\theta \geq \theta_0 | X_1, \ldots, X_n)$

Para a distribuição uniforme no intervalo $[ \max(X_i) - 1, \min(X_i)]$:

Se $\theta_0 < \max(X_i) - 1$, então $ P(_0 | X_1, , X_n) = 1$.

Se $\theta_0 > \min(X_i)$, então $ P(_0 | X_1, , X_n) = 0$.

Se $\max(X_i) - 1 \leq \theta_0 \leq \min(X_i)$:

$P(\theta \geq \theta_0 | X_1, \ldots, X_n) = \frac{\min(X_i) - \theta_0}{\min(X_i) - (\max(X_i) - 1)}$

Questão 2

Considere o modelo $X \sim \text{Uniforme} (0, \theta)$ $(0 < \theta < \infty)$.

Mostre que $\sigma = \theta^{-1}$ é um parâmetro de escala.

Considere $X \sim \text{Uniforme}(0, \theta)$. A função densidade de probabilidade (pdf) de $X$ é:

$f_X(x; \theta) = \theta^{-1} \mathbb{I}_{\{0 \leq x \leq \theta\}}$

Vamos redefinir a distribuição em termos de $\sigma$. Como $\sigma = \theta^{-1}$, então $\theta = \sigma^{-1}$. A pdf em termos de $\sigma$ é:

$f_X(x; \sigma) = \sigma \mathbb{I}_{\{0 \leq x \leq \sigma^{-1}\}}$

Vamos multiplicar $\sigma$ por uma constante $c$. A nova densidade será:

$f_Y(y; c\sigma) = c\sigma \mathbb{I}_{\{0 \leq y \leq (c\sigma)^{-1}\}}$

Essa é uma escala da densidade original, mostrando que $\sigma = \theta^{-1}$ é um parâmetro de escala.

Dada a priori não informativa usual para o modelo de escala, $p(\sigma) \propto \sigma^{-1}$, e uma amostra aleatória $X_1, \ldots, X_n$ do modelo acima, discuta a propriedade da distribuição a posteriori.

Segue que:

$L(\theta | X_1, \ldots, X_n) = \prod_{i=1}^n f_X(x_i; \theta)$

Como $X_i \sim \text{Uniforme}(0, \theta)$, a verossimilhança é:

$L(\theta | X_1, \ldots, X_n) = \theta^{-n} \mathbb{I}_{\{0 \leq \max(X_i) \leq \theta\}}$

A priori não informativa é $p(\sigma) \propto \sigma^{-1}$. Como $\sigma = \theta^{-1}$, temos $p(\theta) \propto \theta^{-1}$. Logo, a distribuição posteriori será:

$p(\theta | X_1, \ldots, X_n) \propto L(\theta | X_1, \ldots, X_n) p(\theta)$

$p(\theta | X_1, \ldots, X_n) \propto \theta^{-n} \mathbb{I}_{\{0 \leq \max(X_i) \leq \theta\}} \theta^{-1}$

$p(\theta | X_1, \ldots, X_n) \propto \theta^{-(n+1)} \mathbb{I}_{\{\theta \geq \max(X_i)\}}$

Por fim:

$p(\theta | X_1, \ldots, X_n) = \frac{(n \max(X_i))^n}{\theta^{n+1}} \mathbb{I}_{\{\theta \geq \max(X_i)\}}$

Esta é uma distribuição de Pareto com forma $n$ e escala $\max(X_i)$.

Na situação da parte (b), ache:
os estimadores bayesianos sob Perda Quadrática e sob Perda Absoluta.

Perda Quadrática: Para a função de perda quadrática $L(\hat{\theta}, \theta) = (\hat{\theta} - \theta)^2$, o estimador bayesiano é a média da distribuição a posteriori. A média de uma distribuição Pareto $\text{Pareto}(\max(X_i), n)$ é dada por: $\hat{\theta}_{\text{PQ}} = \frac{n \max(X_i)}{n - 1}$ , para $n > 1$.
Perda Absoluta: Para a função de perda absoluta $L(\hat{\theta}, \theta) = |\hat{\theta} - \theta|$, o estimador bayesiano é a mediana da distribuição a posteriori. A mediana de uma distribuição Pareto $\text{Pareto}(\max(X_i), n)$ é dada por:

$\hat{\theta}_{\text{PA}} = \max(X_i) \cdot 2^{1/(n)}$

o intervalo HPD para $\theta$ com credibilidade $100(1 - \alpha)\%$.

Para uma distribuição Pareto $\text{Pareto}(\max(X_i), n+1)$, o intervalo HPD $100(1 - \alpha)\%$ pode ser calculado encontrando os quantis apropriados. Seja $F_{\text{Pareto}}^{-1}(p)$ a função quantil da distribuição Pareto:

$\text{HPD} = [ \max(X_i), \max(X_i) (1 - \alpha)^{-\frac{1}{n+1}} ]$

a probabilidade a posteriori do evento $\theta > \theta_0$.

A probabilidade a posteriori do evento $\theta > \theta_0$ é dada pela função de sobrevivência da distribuição Pareto:

$P(\theta > \theta_0 | X_1, \ldots, X_n) = \left( \frac{\max(X_i)}{\theta_0} \right)^{n} \mathbb{I}_{\{\theta_0 \geq \max(X_i)\}}$

Questão 3

Considere o modelo $x_1, \ldots, x_n | \theta \sim \text{Ber}(\theta)$.

Calcule a priori de Jeffreys e mostre que ela é própria.

A priori de Jeffreys é uma escolha não informativa que é proporcional à raiz quadrada do determinante da matriz de informação de Fisher. Para o modelo Bernoulli $\text{Ber}(\theta)$, vamos calcular a matriz de informação de Fisher e encontrar a priori de Jeffreys.

Para uma amostra $x_1, x_2, \ldots, x_n$ de uma distribuição Bernoulli com parâmetro $\theta$, a função de verossimilhança é:

$L(\theta) = \theta^{\sum_{i=1}^n x_i} (1 - \theta)^{n - \sum_{i=1}^n x_i}$

A log-verossimilhança é:

$\log L(\theta) = \left(\sum_{i=1}^n x_i\right) \log \theta + \left(n - \sum_{i=1}^n x_i\right) \log (1 - \theta)$

A informação de Fisher é a expectativa do quadrado da derivada da log-verossimilhança em relação a $\theta$:

$I(\theta) = -\mathbb{E}\left[\frac{\partial^2}{\partial \theta^2} \log L(\theta)\right]$

Calculando a segunda derivada da log-verossimilhança:

$\frac{\partial}{\partial \theta} \log L(\theta) = \frac{\sum_{i=1}^n x_i}{\theta} - \frac{n - \sum_{i=1}^n x_i}{1 - \theta}$

$\frac{\partial^2}{\partial \theta^2} \log L(\theta) = -\frac{\sum_{i=1}^n x_i}{\theta^2} - \frac{n - \sum_{i=1}^n x_i}{(1 - \theta)^2}$

Tomando a expectativa, considerando que $\mathbb{E}[x_i] = \theta$:

$I(\theta) = n \left(\frac{1}{\theta} + \frac{1}{1 - \theta}\right)$

A priori de Jeffreys é proporcional à raiz quadrada da informação de Fisher:

$\pi_J(\theta) \propto \sqrt{I(\theta)}$

$\pi_J(\theta) \propto \sqrt{n \left(\frac{1}{\theta} + \frac{1}{1 - \theta}\right)}$

$\pi_J(\theta) \propto \sqrt{\frac{n}{\theta(1 - \theta)}}$

$\pi_J(\theta) \propto \frac{1}{\sqrt{\theta(1 - \theta)}}$

Esta é a distribuição Beta $\text{Beta}(1/2, 1/2)$, que é própria porque a integral da densidade de Jeffreys sobre o intervalo $[0, 1]$ é finita.

Considere uma amostra aleatória $X_1, \ldots, X_n$ do modelo acima e a priori de Jeffreys. Ache o estimador bayesiano de $\theta$ sob Perda Quadrática e explique como achar um intervalo HPD para $\theta$ com credibilidade $100(1 - \alpha)\%$.

Temos que a verossimilhança é:

$L(\theta | x_1, \ldots, x_n) = \theta^{\sum_{i=1}^n x_i} (1 - \theta)^{n - \sum_{i=1}^n x_i}$

Com a priori de Jeffreys, a posteriori resulta em:

$\pi(\theta | x_1, \ldots, x_n) \propto \theta^{\sum_{i=1}^n x_i} (1 - \theta)^{n - \sum_{i=1}^n x_i} \cdot \frac{1}{\sqrt{\theta(1 - \theta)}}$

$\pi(\theta | x_1, \ldots, x_n) \propto \theta^{\sum_{i=1}^n x_i - 1/2} (1 - \theta)^{n - \sum_{i=1}^n x_i - 1/2}$

Esta é uma distribuição Beta $\text{Beta}(\sum_{i=1}^n x_i + 1/2, n - \sum_{i=1}^n x_i + 1/2)$.

O estimador bayesiano sob perda quadrática é a média da distribuição a posteriori Beta:

$\hat{\theta}_{\text{PQ}} = \mathbb{E}[\theta | x_1, \ldots, x_n] = \frac{\sum_{i=1}^n x_i + 1/2}{n + 1}$

Para encontrar o intervalo HPD, precisamos calcular os quantis da distribuição Beta posteriori $\text{Beta}(\sum_{i=1}^n x_i + 1/2, n - \sum_{i=1}^n x_i + 1/2)$.

Como a questão não explicita o total de observações e a soma das observações, podemos visualizar o comportamento do HPD ao passo em que o tamanho amostral aumenta:

Show/hide code

def hpd_interval(a, b, level=0.99, sample_size=100000):
    sample = beta.rvs(a, b, size=sample_size)
    sample_sorted = np.sort(sample)
    interval_idx_inc = int(np.floor(level * sample_size))
    n_intervals = sample_size - interval_idx_inc
    interval_width = sample_sorted[interval_idx_inc:] - sample_sorted[:n_intervals]
    min_idx = np.argmin(interval_width)
    hpd_min = sample_sorted[min_idx]
    hpd_max = sample_sorted[min_idx + interval_idx_inc]
    return hpd_min, hpd_max

def update_plot(n):
    media_x = 0.7  # Supondo uma média dos xi
    sum_x_i = media_x * n
    a = sum_x_i + 0.5
    b = n - sum_x_i + 0.5
    
    x = np.linspace(0, 1, 1000)
    y = beta.pdf(x, a, b)
    
    hpd_99 = hpd_interval(a, b, level=0.99)
    hpd_99_formatted = (round(hpd_99[0], 3), round(hpd_99[1], 3))
    
    plt.figure(figsize=(6, 3))
    plt.plot(x, y, label=f'Beta({a:.2f}, {b:.2f})')
    plt.fill_between(x, y, where=(x >= hpd_99[0]) & (x <= hpd_99[1]), color='gray', alpha=0.5, label='HPD')
    plt.title(f'Distribuição Beta e Intervalo HPD de 99%\nIntervalo HPD: {hpd_99_formatted}')
    plt.xlabel('θ')
    plt.ylabel('Densidade de Probabilidade')
    plt.legend()
    plt.show()

n_slider = IntSlider(min=1, max=100, step=1, value=10, description='n')

interactive_plot = interactive(update_plot, n=n_slider)

display(interactive_plot)

O gráfico supracitado foi feito pressupondo que a média de $x_i$ fosse 0.7 , pode-se perceber que a média da distribuição posteriori se torna mais estável e menos sensível a pequenas alterações nos dados conforme 𝑛 aumenta. Isso significa que a estimativa pontual de 𝜃 se torna mais confiável.

Os parâmetros $a$ e $b$ da distribuição Beta aumentam linearmente com $n$, tornando a distribuição mais estreita. Para tamanhos amostrais pequenos, a distribuição Beta toma formas largas e até assimétricas, dependendo dos valores observados. Conforme $n$ aumenta, a distribuição tende a se tornar mais simétrica e mais próxima de uma distribuição normal centrada na média posteriori e ,consequentemente, reduzindo o comprimento do intervalo HPD, convergindo para um valor fixo muito próximo da média verdadeira de $\theta$.

Lista 4

Questão 1

Considere uma amostra $y_1, \ldots, y_n$ da distribuição Normal com média $\mu$ e variância $\sigma^2 = 1/\tau$ desconhecidas, e suponha que a priori a distribuição de $(\mu, \tau)$ é especificada da seguinte forma: $\mu | \tau \sim N(\mu_0, 1/\lambda_0 \tau)$ e $\tau \sim \text{Gama}(\alpha_0, \beta_0)$, onde $\lambda_0$, $\alpha_0$ e $\beta_0$ são positivas.

Ache a distribuição a posteriori de $p(\mu, \tau | D)$ e as distribuições marginais $p(\mu | D)$ e $p(\tau | D)$.

Ressaltando:

Priori: $\mu | \tau \sim N(\mu_0, 1/\lambda_0 \tau)$ $\tau \sim \text{Gama}(\alpha_0, \beta_0)$

A função de verossimilhança para uma amostra $y_1, \ldots, y_n$ é:

$p(y | \mu, \tau) = \left( \frac{\tau}{2\pi} \right)^{n/2} \exp\left( -\frac{\tau}{2} \sum_{i=1}^n (y_i - \mu)^2 \right)$

Temos a posteriori:

$p(\mu, \tau | y) \propto p(y | \mu, \tau) p(\mu | \tau) p(\tau)$

Substituindo as expressões:

\[ \begin{align*} p(\mu, \tau \mid y) &\propto \left( \frac{\tau}{2\pi} \right)^{n/2} \exp\left( -\frac{\tau}{2} \sum_{i=1}^n (y_i - \mu)^2 \right) \\ &\quad \times \left( \frac{\lambda_0 \tau}{2\pi} \right)^{1/2} \exp\left( -\frac{\lambda_0 \tau}{2} (\mu - \mu_0)^2 \right) \\ &\quad \times \frac{\beta_0^{\alpha_0}}{\Gamma(\alpha_0)} \tau^{\alpha_0 - 1} \exp(-\beta_0 \tau) \end{align*} \]

Agrupando os termos dependentes de $\mu$ e $\tau$:

$p(\mu, \tau | y) \propto \tau^{(n+1)/2 + \alpha_0 - 1} \exp\left( -\frac{\tau}{2} \left[ \sum_{i=1}^n (y_i - \mu)^2 + \lambda_0 (\mu - \mu_0)^2 \right] - \beta_0 \tau \right)$

Para integrar $\mu$, completamos o quadrado no expoente:

$\sum_{i=1}^n (y_i - \mu)^2 + \lambda_0 (\mu - \mu_0)^2$

Expandindo e agrupando os termos:

$\sum_{i=1}^n y_i^2 - 2\mu \sum_{i=1}^n y_i + n\mu^2 + \lambda_0 \mu_0^2 - 2\lambda_0 \mu \mu_0 + \lambda_0 \mu^2$

Agrupando os termos quadráticos, lineares e constantes em $\mu$:

$\left( n + \lambda_0 \right) \mu^2 - 2 \left( \sum_{i=1}^n y_i + \lambda_0 \mu_0 \right) \mu + \sum_{i=1}^n y_i^2 + \lambda_0 \mu_0^2$

Completando o quadrado:

$\left( n + \lambda_0 \right) \left( \mu - \frac{\sum_{i=1}^n y_i + \lambda_0 \mu_0}{n + \lambda_0} \right)^2 + \left( \sum_{i=1}^n y_i^2 + \lambda_0 \mu_0^2 \right) - \frac{\left( \sum_{i=1}^n y_i + \lambda_0 \mu_0 \right)^2}{n + \lambda_0}$

Assim, o expoente se torna:

$-\frac{\tau}{2} \left[ (n + \lambda_0) \left( \mu - \frac{\sum_{i=1}^n y_i + \lambda_0 \mu_0}{n + \lambda_0} \right)^2 + \sum_{i=1}^n y_i^2 + \lambda_0 \mu_0^2 - \frac{\left( \sum_{i=1}^n y_i + \lambda_0 \mu_0 \right)^2}{n + \lambda_0} \right]$

A integral em $\mu$ da expressão acima é uma integral gaussiana, logo:

$\int_{-\infty}^{\infty} \exp \left( -\frac{\tau (n + \lambda_0)}{2} \left( \mu - \frac{\sum_{i=1}^n y_i + \lambda_0 \mu_0}{n + \lambda_0} \right)^2 \right) d\mu = \sqrt{\frac{2\pi}{\tau (n + \lambda_0)}}$

Portanto, ao desprezarmos essa constante, a distribuição marginal de $\tau$ é:

$p(\tau | y) \propto \tau^{(n+1)/2 + \alpha_0 - 1} \exp \left( -\tau \left( \frac{\sum_{i=1}^n y_i^2 + \lambda_0 \mu_0^2 - \frac{\left( \sum_{i=1}^n y_i + \lambda_0 \mu_0 \right)^2}{n + \lambda_0}}{2} + \beta_0 \right) \right)$

Simplificando a constante no expoente:

$p(\tau | y) \propto \tau^{(n+1)/2 + \alpha_0 - 1} \exp \left( -\tau \left( \frac{\sum_{i=1}^n (y_i - \bar{y})^2}{2} + \frac{\lambda_0 n (\bar{y} - \mu_0)^2}{2(\lambda_0 + n)} + \beta_0 \right) \right)$

Reconhecendo a forma da distribuição Gama, temos:

$\tau | y \sim \text{Gama} \left( \alpha_0 + \frac{n}{2}, \beta_0 + \frac{1}{2} \left[ \sum_{i=1}^n (y_i - \bar{y})^2 + \frac{\lambda_0 n (\bar{y} - \mu_0)^2}{\lambda_0 + n} \right] \right)$

Dada a $\tau$, a distribuição condicional de $\mu$ é:

$\mu | \tau, y \sim N \left( \frac{\lambda_0 \mu_0 + n \bar{y}}{\lambda_0 + n}, \frac{1}{\tau (\lambda_0 + n)} \right)$

$\tau | y \sim \text{Gama} \left( \alpha_0 + \frac{n}{2}, \beta_0 + \frac{1}{2} \left[ \sum_{i=1}^n (y_i - \bar{y})^2 + \frac{\lambda_0 n (\bar{y} - \mu_0)^2}{\lambda_0 + n} \right] \right)$

Dada a $\tau$, a distribuição condicional de $\mu$ é:

$\mu | \tau, y \sim N \left( \frac{\lambda_0 \mu_0 + n \bar{y}}{\lambda_0 + n}, \frac{1}{\tau (\lambda_0 + n)} \right)$

Discuta neste contexto o uso da distribuição a priori não informativa $p(\mu, \tau) \propto 1/\tau$. A distribuição a posteriori é própria? Qual é a relação da distribuição $p(\mu | D)$ e os resultados clássicos?

Essa priori é escolhida para ser não informativa, significando que temos pouca ou nenhuma informação a priori sobre $\mu$ e $\tau$.

Distribuição a posteriori: Com a priori não informativa, a posteriori ainda é própria, pois as integrais de normalização convergem.
Comparação com Resultados Clássicos:
- A distribuição marginal de $\mu$ se aproxima da distribuição normal clássica para a média de uma amostra, especialmente à medida que $n$ aumenta.
- A distribuição marginal de $\tau$ se aproxima da distribuição qui-quadrado clássica (quando reescalada), que é usada para inferência sobre a variância na estatística frequentista.

Questão 2c

Seja $n = (n_1, \ldots, n_k)$ um vetor aleatório com distribuição multinomial e densidade $p(n | \theta) \propto \prod_{i=1}^k \theta_i^{n_i}$, onde $\theta = (\theta_1, \ldots, \theta_k)$, $\theta_i > 0$ e $\sum_{i=1}^k \theta_i = 1$. Considere a priori para $\theta$ uma distribuição de Dirichlet com parâmetro $\alpha = (\alpha_1, \ldots, \alpha_k)$, isto é $p(\theta) \propto \prod_{i=1}^k \theta_i^{\alpha_i - 1}$.

No caso particular $\alpha = (0, 0, \ldots, 0)$ a distribuição a priori de $\theta$ é imprópria. Mostre que a distribuição a posteriori é própria se e somente se $n_i > 0$ para $i = 1, 2, \ldots, k$.

Quando $\alpha = (0, 0, \ldots, 0)$, a distribuição a priori de $\theta$ é:

$p(\theta) \propto \prod_{i=1}^k \theta_i^{-1}$

Esta distribuição é imprópria porque a integral sobre o simplex $\{\theta \in \mathbb{R}^k : \theta_i > 0, \sum_{i=1}^k \theta_i = 1\}$ não converge.

A função de verossimilhança da distribuição multinomial é:

$p(n | \theta) \propto \prod_{i=1}^k \theta_i^{n_i}$

Temos assim a posteriori:

$p(\theta | n) \propto p(n | \theta) p(\theta) \propto \left( \prod_{i=1}^k \theta_i^{n_i} \right) \left( \prod_{i=1}^k \theta_i^{-1} \right) = \prod_{i=1}^k \theta_i^{n_i - 1}$

Esta é a forma da distribuição Dirichlet com parâmetros $(n_1, n_2, \ldots, n_k)$:

$p(\theta | n) \propto \prod_{i=1}^k \theta_i^{n_i - 1}$

Para que a distribuição a posteriori seja própria, a integral da densidade da posteriori sobre o simplex deve ser finita. A distribuição Dirichlet $\text{Dirichlet}(n_1, n_2, \ldots, n_k)$ é própria se e somente se todos os $n_i > 0$. Se qualquer $n_i = 0$, o termo $\theta_i^{-1}$ aparece na densidade, e a integral não é finita.

Portanto, a distribuição a posteriori é própria se e somente se $n_i > 0$ para todos os $i = 1, 2, \ldots, k$.

Questão 3

Na véspera do primeiro turno para a eleição de governador do DF de 2010, a Datafolha divulgou uma pesquisa indicando que, de 891 eleitores entrevistados que já tinham decidido em quem votar, Agnelo Queiroz tinha a preferência de 467, Weslian Roriz a de 315 e outros candidatos a de 109 eleitores.

Formule um modelo para analisar esses dados. O interesse centra fundamentalmente em três perguntas:

A eleição poderia ser definida no primeiro turno?

Usaremos um modelo baseado na distribuição Multinomial para representar os votos dos eleitores. A distribuição Multinomial é adequada aqui porque modela a probabilidade de um certo número de sucessos em várias categorias.

Seja $\theta = (\theta_1, \theta_2, \theta_3)$ as proporções de votos para Agnelo, Weslian e outros candidatos, respectivamente, onde:

$\theta_1 + \theta_2 + \theta_3 = 1$

Os dados $n = (n_1, n_2, n_3)$ seguem uma distribuição multinomial:

$(n_1, n_2, n_3) \sim \text{Multinomial}(n, \theta)$

Para a análise bayesiana, assumimos uma distribuição a priori para $\theta$. Podemos usar uma distribuição de Dirichlet com parâmetros $\alpha = (\alpha_1, \alpha_2, \alpha_3)$:

$\theta \sim \text{Dirichlet}(\alpha)$

Dado que os dados seguem uma distribuição multinomial e a priori é uma Dirichlet, a distribuição a posteriori de $\theta$ também é uma Dirichlet com parâmetros $\alpha' = (\alpha_1 + n_1, \alpha_2 + n_2, \alpha_3 + n_3)$:

$\theta | n \sim \text{Dirichlet}(\alpha_1 + n_1, \alpha_2 + n_2, \alpha_3 + n_3)$

Para uma priori não informativa, podemos assumir $\alpha_i = 1$ para todos os $i$.

Para que a eleição seja definida no primeiro turno, um candidato deve obter mais de 50% dos votos válidos. Vamos calcular a probabilidade posteriori de $\theta_1 > 0.5$:

$P(\theta_1 > 0.5 | n)$

O candidato Agnelo poderia ser eleito no primeiro turno?

Semelhante à (a), precisamos calcular:

$P(\theta_1 > 0.5 | n)$

Qual será a diferença na porcentagem dos votos válidos entre os dois primeiros colocados?

Calculamos a diferença na proporção dos votos válidos entre Agnelo e Weslian:

$P(|\theta_1 - \theta_2| > d | n)$

onde $d$ é a diferença que desejamos calcular.

Podemos implementar isso usando simulações a partir da distribuição posteriori Dirichlet para estimar essas probabilidades.

Show/hide code

n_total = 891
n_agnelo = 467
n_weslian = 315
n_outros = 109

alpha = np.array([1, 1, 1])

alpha_post = alpha + np.array([n_agnelo, n_weslian, n_outros])

n_samples = 100000
samples = dirichlet(alpha_post).rvs(n_samples)

prob_agnelo_eleito = np.mean(samples[:, 0] > 0.5)
dif_agnelo_weslian = samples[:, 0] - samples[:, 1]
diff_percent = np.mean(dif_agnelo_weslian)

print(f"Probabilidade de Agnelo ser eleito no primeiro turno: {prob_agnelo_eleito:.2f}")
print(f"Diferença média na porcentagem entre Agnelo e Weslian: {diff_percent:.2f}")

plt.figure(figsize=(5,5))
x = np.linspace(0, 1, 1000)
plt.plot(x, beta.pdf(x, alpha_post[0], alpha_post[1] + alpha_post[2]), label="Agnelo")
plt.plot(x, beta.pdf(x, alpha_post[1], alpha_post[0] + alpha_post[2]), label="Weslian")
plt.axvline(0.5, color='red', linestyle='--', label="50% Threshold")
plt.xlabel("Proporção")
plt.ylabel("Densidade de Probabilidade")
plt.title("Distribuições Posteriori das Proporções de Votos")
plt.legend(bbox_to_anchor=(1.05, 1), loc='upper right')
plt.grid(True)
plt.tight_layout() 
plt.show()

Probabilidade de Agnelo ser eleito no primeiro turno: 0.92
Diferença média na porcentagem entre Agnelo e Weslian: 0.17

Com base nos dados da pesquisa e na análise bayesiana, é altamente provável que Agnelo Queiroz seja eleito no primeiro turno, com uma vantagem média de 17% sobre Weslian Roriz. A probabilidade de Agnelo obter mais de 50% dos votos válidos é muito alta (91.97%), o que sugere que a eleição pode ser decidida no primeiro turno.

Lista 5

Questão 2

Considere $n = 12$ ensaios de Bernoulli com probabilidade de sucesso $p$ e $y = 9$ sucessos. Suponha que a distribuição a priori de $p$ é especificada de forma que $\eta = \log \left( \frac{p}{1 - p} \right)$ segue uma distribuição Normal com média $\mu = 0$ e variância $\sigma^2 = 1$.

De imediato temos:

$\eta \sim N(0, \sigma^2)$

O parâmetro $p$ pode ser expresso em termos de $\eta$ como:

$p = \frac{e^\eta}{1 + e^\eta}$

A verossimilhança para a distribuição Bernoulli é:

$p(y | p) = \binom{n}{y} p^y (1 - p)^{n - y}$

Para $n = 12$ e $y = 9$:

$p(9 | p) = \binom{12}{9} p^9 (1 - p)^3$

Para calcular $\mathbb{E}(p | y = 9)$ e $\text{DP}(p | y = 9)$, podemos usar a aproximação de Laplace, que envolve encontrar a moda da posteriori e aproximar a distribuição posteriori por uma Normal em torno da moda.

Podemos implementar isso em Python para calcular as aproximações solicitadas.

Show/hide code

n = 12
y = 9
mu = 0

variances = [1, 4, 9]

def posterior_approximation(y, n, mu, sigma2):

    # Log-verossimilhança em termos de eta
    def log_likelihood_eta(eta):
        p = expit(eta)  
        return y * np.log(p) + (n - y) * np.log(1 - p)
    
    def prior_eta(eta):
        return norm.logpdf(eta, loc=mu, scale=np.sqrt(sigma2))
    
    def log_posterior_eta(eta):
        return log_likelihood_eta(eta) + prior_eta(eta)
    
    def dlog_posterior_eta(eta):
        p = expit(eta)
        return y - n * p - eta / sigma2
    
    def d2log_posterior_eta(eta):
        p = expit(eta)
        return -n * p * (1 - p) - 1 / sigma2
    
    # Modo da posteriori (máximo a posteriori)
    from scipy.optimize import minimize
    result = minimize(lambda eta: -log_posterior_eta(eta), x0=0, method='Newton-CG', jac=lambda eta: -dlog_posterior_eta(eta), hess=lambda eta: -d2log_posterior_eta(eta))
    eta_mode = result.x[0]
    
    # Variância em torno da moda
    eta_var = -1 / d2log_posterior_eta(eta_mode)
    
    eta_samples = np.random.normal(loc=eta_mode, scale=np.sqrt(eta_var), size=10000)
    p_samples = expit(eta_samples)
    
    mean_p = np.mean(p_samples)
    std_p = np.std(p_samples)
    prob_p_greater_half = np.mean(p_samples > 0.5)
    
    return p_samples, mean_p, std_p, prob_p_greater_half

fig, axes = plt.subplots(len(variances), 1, figsize=(5,10))

for ax, sigma2 in zip(axes, variances):
    p_samples, mean_p, std_p, prob_p_greater_half = posterior_approximation(y, n, mu, sigma2)
    
    ax.hist(p_samples, bins=50, density=True, alpha=0.75)
    ax.axvline(mean_p, color='r', linestyle='--', linewidth=2, label=f'Média = {mean_p:.3f}')
    ax.axvline(0.5, color='g', linestyle='-.', linewidth=2, label='$p = 0.5$')
    ax.set_title(f'Posterior de $p$ para $\sigma^2 = {sigma2}$', fontsize=14)
    ax.set_xlabel('$p$', fontsize=12)
    ax.set_ylabel('Densidade', fontsize=12)
    ax.legend(loc='upper left', fontsize=12)
    ax.grid(True)
    
    
    print(f"Para σ² = {sigma2}:")
    print(f"  E(p | y = 9) = {mean_p:.4f}")
    print(f"  DP(p | y = 9) = {std_p:.4f}")
    print(f"  Pr(p > 0.5 | y = 9) = {prob_p_greater_half:.4f}")
    print()

plt.tight_layout()
plt.show()

Para σ² = 1:
  E(p | y = 9) = 0.6736
  DP(p | y = 9) = 0.1099
  Pr(p > 0.5 | y = 9) = 0.9275

Para σ² = 4:
  E(p | y = 9) = 0.7144
  DP(p | y = 9) = 0.1168
  Pr(p > 0.5 | y = 9) = 0.9502

Para σ² = 9:
  E(p | y = 9) = 0.7240
  DP(p | y = 9) = 0.1208
  Pr(p > 0.5 | y = 9) = 0.9469

Percebemos que a medida que a variância da priori aumenta, a incerteza na posteriori também aumenta, refletindo uma maior incerteza a priori sobre $p$, sendo propagada para a distribuição a posteriori.